home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / regress1.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  120 lines
  1. ; $Id: regress1.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1993-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6. FUNCTION REGRESS1,X,Y,W,YFIT,A0,SIGMA,FTEST,R,RMUL,CHISQ
  7.  
  8. ;+
  9. ; NAME:
  10. ;    REGRESS1
  11. ;
  12. ; PURPOSE:
  13. ;    Multiple linear regression fit.
  14. ;    Fit the function:
  15. ;    Y(i) = A0 + A(0)*X(0,i) + A(1)*X(1,i) + ... + 
  16. ;        A(Nterms-1)*X(Nterms-1,i)
  17. ;
  18. ; CATEGORY:
  19. ;    G2 - Correlation and regression analysis.
  20. ;
  21. ; CALLING SEQUENCE:
  22. ;    Result = REGRESS(X, Y, W, YFIT, A0, SIGMA, FTEST, R, RMUL, CHISQ)
  23. ;
  24. ; INPUTS:
  25. ;    X:    array of independent variable data.  X must 
  26. ;        be dimensioned (Nterms, Npoints) where there are Nterms 
  27. ;        coefficients to be found (independent variables) and 
  28. ;        Npoints of samples.
  29. ;
  30. ;    Y:    vector of dependent variable points, must have Npoints 
  31. ;        elements.
  32. ;
  33. ;       W:    vector of weights for each equation, must be a Npoints 
  34. ;        elements vector.  For instrumental weighting 
  35. ;        w(i) = 1/standard_deviation(Y(i)), for statistical 
  36. ;        weighting w(i) = 1./Y(i).   For no weighting set w(i)=1,
  37. ;        and also set the RELATIVE_WEIGHT keyword.
  38. ;
  39. ; OUTPUTS:
  40. ;    Function result = coefficients = vector of 
  41. ;    Nterms elements.  Returned as a column vector.
  42. ;
  43. ; OPTIONAL OUTPUT PARAMETERS:
  44. ;    Yfit:    array of calculated values of Y, Npoints elements.
  45. ;
  46. ;    A0:    Constant term.
  47. ;
  48. ;    Sigma:    Vector of standard deviations for coefficients.
  49. ;
  50. ;    Ftest:    value of F for test of fit.
  51. ;
  52. ;    Rmul:    multiple linear correlation coefficient.
  53. ;
  54. ;    R:    Vector of linear correlation coefficient.
  55. ;
  56. ;    Chisq:    Reduced weighted chi squared.
  57. ;
  58. ; KEYWORDS:
  59. ;RELATIVE_WEIGHT: if this keyword is non-zero, the input weights
  60. ;        (W vector) are assumed to be relative values, and not based
  61. ;        on known uncertainties in the Y vector.    This is the case for
  62. ;        no weighting W(*) = 1.
  63. ;
  64. ; PROCEDURE:
  65. ;    Adapted from the program REGRES, Page 172, Bevington, Data Reduction
  66. ;    and Error Analysis for the Physical Sciences, 1969.
  67. ;
  68. ; MODIFICATION HISTORY:
  69. ;    Written, DMS, RSI, September, 1982.
  70. ;    Added RELATIVE_WEIGHT keyword, W. Landsman, August 1991
  71. ;-
  72. ;
  73. On_error,2              ;Return to caller if an error occurs 
  74. SY = SIZE(Y)            ;Get dimensions of x and y.  
  75. SX = SIZE(X)
  76. IF (N_ELEMENTS(W) NE SY(1)) OR (SX(0) NE 2) OR (SY(1) NE SX(2)) THEN $
  77.   message, 'Incompatible arrays.'
  78. ;
  79. NTERM = SX(1)           ;# OF TERMS
  80. NPTS = SY(1)            ;# OF OBSERVATIONS
  81. ;
  82. SW = TOTAL(W)           ;SUM OF WEIGHTS
  83. YMEAN = TOTAL(Y*W)/SW   ;Y MEAN
  84. XMEAN = (X * (REPLICATE(1.,NTERM) # W)) # REPLICATE(1./SW,NPTS)
  85. WMEAN = SW/NPTS
  86. WW = W/WMEAN
  87. ;
  88. NFREE = NPTS-1          ;DEGS OF FREEDOM
  89. SIGMAY = SQRT(TOTAL(WW * (Y-YMEAN)^2)/NFREE) ;W*(Y(I)-YMEAN)
  90. XX = X- XMEAN # REPLICATE(1.,NPTS)      ;X(J,I) - XMEAN(I)
  91. WX = REPLICATE(1.,NTERM) # WW * XX      ;W(I)*(X(J,I)-XMEAN(I))
  92. SIGMAX = SQRT( XX*WX # REPLICATE(1./NFREE,NPTS)) ;W(I)*(X(J,I)-XM)*(X(K,I)-XM)
  93. R = WX #(Y - YMEAN) / (SIGMAX * SIGMAY * NFREE)
  94.  
  95.  
  96. WW1 = WX # TRANSPOSE(XX)
  97. d = determ(WW1)
  98. if d eq 0 THEN return, 1.e+30
  99. if d lt 1.e-13 THEN BEGIN 
  100.   print,"regression-- Determinant of correlation matrix less than 1.e-13."
  101.   print,"             Hit control C to terminate."
  102. ENDIF
  103.  
  104. ARRAY = INVERT((WX # TRANSPOSE(XX))/(NFREE * SIGMAX #SIGMAX))
  105. A = (R # ARRAY)*(SIGMAY/SIGMAX)         ;GET COEFFICIENTS
  106. YFIT = A # X                            ;COMPUTE FIT
  107. A0 = YMEAN - TOTAL(A*XMEAN)             ;CONSTANT TERM
  108. YFIT = YFIT + A0                        ;ADD IT IN
  109. FREEN = NPTS-NTERM-1 > 1                ;DEGS OF FREEDOM, AT LEAST 1.
  110.  
  111. CHISQ = TOTAL(WW*(Y-YFIT)^2)*WMEAN/FREEN ;WEIGHTED CHI SQUARED
  112. IF KEYWORD_SET(relative_weight) then varnce = chisq $
  113.                                 else varnce = 1./wmean
  114. sigma = sqrt(array(indgen(nterm)*(nterm+1))*varnce/(nfree*sigmax^2)) ;Error term
  115. RMUL = TOTAL(A*R*SIGMAX/SIGMAY)         ;MULTIPLE LIN REG COEFF
  116. IF RMUL LT 1. THEN FTEST = RMUL/NTERM / ((1.-RMUL)/FREEN) ELSE FTEST=1.E6
  117. RMUL = SQRT(RMUL)
  118. RETURN,A
  119. END
  120.